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Abstract — The problem of distributed estimation of a paramet- 
ric physical field is stated as a maximum likelihood estimation 
problem. Sensor observations are distorted by additive white 
Gaussian noise. Prior to data transmission, each sensor quan- 
tizes its observation to M levels. The quantized data are then 
communicated over parallel additive white Gaussian channels to 
a fusion center for a joint estimation. An iterative expectation- 
maximization (EM) algorithm to estimate the unknown param- 
eter is formulated, and its linearized version is adopted for 
numerical analysis. The numerical examples are provided for the 
case of the field modeled as a Gaussian bell. The dependence of 
the integrated mean-square error on the number of quantization 
levels, the number of sensors in the network and the SNR in 
observation and transmission channels is analyzed. 

Index Terms — Distributed estimation, expectation- 
maximization algorithm, maximum likelihood estimation, 
distributed sensor network, sparse data 

I. Introduction 

Distributed sensor networks provide a platform for many 
military and civilian applications. Examples include surveil- 
lance, monitoring wildlife, or controlling the power grid. 
Sensor networks built for these applications are intended to 
solve various problems such as detecting, tracking, classifying, 
counting, and estimating. They are also required to adhere to 
a number of physical constraints including power, bandwidth, 
latency, and complexity. Much research has been reported 
on each of these topics over the past two decades. In the 
field of distributed estimation, for example, various estimation 
problems have been formulated and solved. Many works 
choose either to optimize a distributed sensor network with 
respect to energy consumption during transmission fTlL 151 or 
impose bandwidth constraints and thus focus on designing an 
optimal quantization strategy for the distributed network 0, 
l4l . There are few that involve both constraints (see as an 
example). 

Among research groups working on the problem of dis- 
tributed estimation, there are a few dealing with distributed 
estimation of a field (a multidimensional function, in gen- 
eral) El, 0, 0. Since in many real- world applications 
distributed estimation of a multidimensional function may 
provide additional information that aids in making a high- 
fidelity decision or in solving another inference problem, we 
contribute to this topic by formulating and solving the problem 

N. A. Schmid, M. Alkhweldi, and M. C. Valenti are with the Department 
of Computer Science and Electrical Engineering, West Virginia Univer- 
sity, Morgantown, WV, 26506 USA e-mail: Natalia.Schmid@mail.wvu.edu, 
malkhwel@mix.wvu.edu, and Matthew.Valenti@mail.wvu.edu. 

This work was sponsored by the Office of Naval Research under Award 
No. N00014-09-1-1189. 



of a parametric field estimation from sparse noisy sensor 
measurements. Distributed target localization is another active 
area of research [4 ]. An iterative solution to this problem is a 
second contribution of our work. 

In this paper, the problem of distributed estimation of a 
physical field from sensory data collected by a homogeneous 
sensor network is stated as a maximum likelihood estimation 
problem. The physical field is a deterministic function and has 
a known spatial distribution parameterized by a set of unknown 
parameters, such as the location of an object generating the 
field and the strength of the field in the region occupied 
by the sensors. Sensor observations are distorted by additive 
white Gaussian noise. Prior to transmission, each sensor 
quantizes its observation to M levels. The quantized data 
are then communicated over parallel additive white Gaussian 
channels to the fusion center where the unknown parameters 
of the underlying physical field are estimated. An iterative 
expectation-maximization (EM) algorithm to estimate the un- 
known parameter is formulated, and a simplified numerical 
solution involving additional approximations is developed. The 
numerical examples illustrate the developed approach for the 
case of the field modeled as a Gaussian bell. 

The remainder of the paper is organized as follows. Sec. [TT] 
formulates the problem. Sec. [HT| develops an EM solution. Sec. 
[IV] provides numerical performance evaluation. The summary 
of the developed results is provided in Sec. (Vj 

II. Problem Statement 

Consider a distributed network of homogeneous sensors 
monitoring the environment for the presence of a substance 
or an object. Assume that each substance or object is charac- 
terized by a location parameter and by a spatially distributed 
physical field generated by it. As an example, a ferromagnetic 
object can be viewed as a single or a collection of dipoles 
characterized by a magnetic field that they generate. This 
field can be sensed by a network of magnetometers placed 
in the vicinity of the object. Depending on the design of the 
magnetometers, they may take measurements of a directional 
complex valued magnetic field or of the magnitude of the field 
only. The field generated by a dipole decays as a function 
of the inverse cube of the distance to the dipole. The sensor 
network does not know a priori the location of the dipole as 
well as the type and size of the object. However, the type and 
size of an object can be associated with the strength of the 
magnetic field. Examples of other physical fields include (1) a 
radioactive field that can be modeled as a stationary spatially 
distributed Poisson field with a two-dimensional intensity 
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Fig. 1. Block-diagram of the distributed sensor network. 



function decaying according to the inverse- square law or (2) a 
distribution of pollution or chemical fumes that, if stationary, 
can often be modeled as a Gaussian bell. 

Consider a network of K sensors distributed over an area 
A. The network is calibrated in the sense that the relative 
locations of the sensors are known. Sensors act independently 
of one another and take noisy measurements of a physical 
field G{x,y). A sample of G{x,y) at a location (xk,yk) is 
denoted as Gk = G(xk,yk)- The parametric field G(x,y) is 
characterized by L unknown parameters 6 = . . . , 9l] t . 
To emphasize this dependence we will use both Gk and 
G(xk,Vk : 0) throughout the text. The sensor noise, de- 
note it by Wk : k = 1,...,K is known and modeled as 
Gaussian distributed with mean zero and variance a 2 . The 
noise of sensors is independent and identically distributed. 
Let Rk, k = 1, . . . , K be the noisy samples of the field at 
the location of distributed sensors. Then Rk is modeled as 
R k = G(x k ,y k ) + W k . 

Due to constraints that are imposed by practical technology, 
each sensor may be required to quantize its measurements 
prior to transmitting them to the fusion center (FC). Assume 
that a deterministic quantizer with M quantization levels is 
involved. Let v\ , V2 , • • • , vm be known reproduction points of 
the quantizer. Denote by q(Rk) = qk the quantized version of 
the measurement by the k-th sensor. These data are modulated 
using a digital modulation scheme and then transmitted to the 
FC over noisy parallel channels. The noise in channels is due 
to quantization error and channel impairments, denote it by 
iVfc, k = 1,...,K. Denote by ra(-) a modulation function 
and by d(-) a demodulation function. Let Zi, . . . , Zk, be 
noisy observations received by the FC. Then each Z k is given 
by Zk = d(m{qk)) + iVfc, k = 1, . . . , K. In this work we 
assume that ra(-) and d(-) are linear and that the demodulator 
recovers the quantized signal by using a soft thresholding 
rule. These assumptions allow Zk be approximated by its 
asymptotic counterpart Z k = qk + where N k is a white 
Gaussian noise with variance rf . 

Given the noisy measurements and the relative location of 
the sensors in the network, the task of the FC is to estimate 



the vector parameter 6. A block diagram of the distributed 
network used for estimation of parameters of a physical field 
is shown in Fig. [T] 

In this work we adopt a maximum likelihood (ML) estima- 
tion approach to solve the problem of distributed parameter 
estimation. The joint likelihood function of the independent 
quantized noisy measurements Z\ , Z2 , • • • , Zk can be written 
as 
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where Z is the vector of measurements [Zi, Z2, . . . , Zk} 7 ', 
Pkj are the probabilities for the output of the sensor k to be 
mapped to the j-th reproduction point during the encoding 
process 
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Tj and Tj+i, j = 1,...,M are the boundaries of the j- 
th quantization region. The ML solution is the solution 
that maximizes the expression ([T]). For a numerical example 
in Sec. [TV) the field is modeled as a Gaussian bell with 
three unknown parameters: the strength of the field fi and the 
location parameter (x c ,y c ). 

III. Iterative Solution 

Since the expression for the log-likelihood function ([T]) 
is highly nonlinear in unknown parameters, we develop an 
iterative solution to the problem. We first formulate a set 
of Expectation-Maximization (EM) iterations |9] and then 
involve a Newton's linearization to solve for the unknown 
parameters. 

A. Expectation Maximization Solution 

We select the pairs of random variables (Rk,Nk), k — 
1,2, as complete data. The complete data log- 

likelihood, denote it by l c d(')> is given by 
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The measurements Zi, i = 1 , . . . , K, form incomplete data. 
The mapping from complete data space to incomplete data 
space is given by Zk = q(Rk) + Nk, where q(.) is a known 
quantization function. 

Denote by 0^ an estimate of the vector obtained at 
the /c-th iteration. To update the estimates of parameters we 
alternate the expectation and maximization steps. During the 
expectation step, we evaluate the conditional expectation of 
the complete data log-likelihood: 
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where the expectation is with respect to the conditional 
probability density function of the complete data, given the 
incomplete data (measurements) and the estimates of the 
parameters at the k-th iteration. During the maximization step 
we maximize ((3): 
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To find the conditional expectation we note that the condi- 
tional probability density function (p.d.f.) of Z iy given R i: is 
Gaussian with mean q(Ri) and variance rf and the p.d.f. of 
Ri is Gaussian with mean Gi and variance a 2 . We also note 
that at the k-th iteration the conditional pdf of Ri, given 
implicitly involves the estimates of the parameters obtained at 
the /c-th iteration. 

Denote by the estimate of the field G(x,y) at the 
location (xi,yi) with the vector of parameters 8 replaced by 
their estimates 8^ k \ Then the final expression for the iterative 
evaluation of the unknown parameters can be written as 
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The expression Q(-) is used to denote the Q-function. 

5. Linearization 

The equations ([5]) are nonlinear in and have to be 

solved numerically for each iteration. To simplify the solution, 
we linearize the expression in §5§ by means of Newton's 
method. Denote by F(#( fe+1 )) the vector form of the left 
side in ([5}, which is a mapping from 8^ k+1 ^ to the range of 
F((9^ +1 )). Let J (0 { n +1) ^j be the Jacobian of the mapping. 
The index n indicates the iteration of the Newton's solution. 
Then 8^ solves the following linearized equation: 
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Fig. 2. Squared Gaussian field located at (x c ,y c ) = (4,4) with peak 
parameter 4 and variance 4. 



IV. Numerical Analysis 

In this section, the performance of the distributed ML 
estimator in §5§ is demonstrated on simulated data. A dis- 
tributed network of K sensors is formed by positioning sensors 
at random over an area A of size 8x8. The location of 
each sensor is noted. A Gaussian field shown in Fig. [2] is 
sampled at the location of the i-th sensor, i = 1, . . . , if and 
a sample of randomly generated Gaussian noise with mean 
zero and variance a 2 is added to each field measurement. In 
our simulations, K is varied from 5 to 200 and a 2 is selected 
such that the total signal-to-noise ratio (SNR) of the local 
observations defined as 
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is 15 dB. Each sensor observation is quantized to one of 
M levels using a uniform deterministic quantizer. We set the 
number of quantization levels to M = 8 and the quantization 
step to 8. if parallel white Gaussian noise channels add 
samples of noise with variance r] 2 selected to set the total 
SNR during data transmission defined as 
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to 15 dB, and the FC observes the noisy quantized samples 
of the field. The function q(R(x, y : 8)) in ( 10) is a quantized 
version of R(x, y : 8) = G(x, y : 8) +W. 

First, we illustrate the convergence of the EM algorithm. 
The value of the ML estimate as a function of iteration is 
shown in Figs. |4| [5] and [6] for the peak value of the field, 
its ^-location and its ^-location, respectively. This illustration 
is based on a single realization of the distributed network 
with K = 10 and M = 8. We can observe that with the 
initial values 9 for the peak of the field, 3 for the x-location 
and 3 for the y-location, the algorithm takes about 600 EM 
iterations to converge to the final values 7.90, 3.88, and 3.88, 
respectively. The true values of these parameters are 8, 4, and 
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Fig. 3. The squared difference between the original and reconstructed fields. 
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Fig. 6. Illustration of the EM convergence for M ■ 



Y-location. 




NUMBER OF SENSORS 

Fig. 7. A box plot of the SE between the estimated and true location of the 
object displayed as a function of the number of sensors distributed over the 
area A. The number of quantization levels is set to M = 8. 



Fig. 4. Illustration of the EM convergence for M = 8. Peak parameter. 
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Fig. 5. Illustration of the EM convergence for M = 8. X-location. 



4. The discrepancy between the vectors of estimated and true 
parameters are due to a low sensor density in the network, 
relatively rough quantization, and the distortions due to sensor 
and channel noise. 

The square distance per pixel between the original and 
reconstructed Gaussian fields is displayed in Fig. [3] 

To further analyze the estimation performance, we evaluate 
the mean square error (MSE) between the estimated and true 
location parameters. The MSE is evaluated numerically by 
means of 1000 Monte Carlo simulations. Each vector of esti- 
mated parameters is substituted back in the expression for the 
parametric field, and an integrated square error (ISE) between 
the true and estimated fields is evaluated. The integrated square 
error (ISE) is defined as 
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Fig. 8. Dependence of the simulated ISE on the number of sensors distributed 
over the area A for the case of M = 8. 
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The ISE statistically averaged over 1000 Monte Carlo simu- 
lations is an approximation to integrated mean square error 
(IMSE). 

The dependence of the MSE on the number of sensors, K, 
in the distributed network for the case of M = 8 quantization 
levels is shown in Fig. [7] The dependence of the IMSE on the 
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Fig. 9. Probability of outliers P utUers( T ) — P[SE > r] (expressed 
in percents) as a function of r. The plot is based on 1000 Monte Carlo 
simulations. 
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Fig. 11. A box plot of the SE between the estimated and true location of 
the object displayed as a function of the number of sensors distributed over 
the area A. The number of quantization levels is set to M = 16. 
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Fig. 10. Probability of outliers (expressed in percents) as a function of 
the threshold for different values of SNR in observation and transmission 
channels. 



number of sensors (sensor density) in the distributed network 
for the same value of M is displayed in Fig. [8] The number of 
sensors distributed over the area A is varied from 5 to 200 with 
the step 5. Each box in Fig. [7] and Fig. [8] is generated using 
1000 Monte Carlo realizations of the network and EM runs. 
The central mark in each box is the median. The edges of the 
box present the 25th and 75th percentiles. The dashed vertical 
lines mark the data that extend beyond the two percentiles, but 
not considered as outliers. The outliers are plotted individually 
and marked with a "+" sign. The percentage of outliers due 
to divergence of the EM algorithm is depicted in Fig. [9] Note 
the large percentage of outliers for small values of K, K = 
10, 15, 20. These correspond to the case when one of the three 
parameters did not converge to its true value. 

The results indicate that the location estimation and the field 
reconstruction of a relatively good quality is possible with 
M = 8 and the number of sensors equal or exceeding 20. 

Fig. [TO] compares the percentage of outliers plotted as a 
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Fig. 12. Dependence of the simulated ISE on the number of sensors 
distributed over the area A. The number of quantization levels is set to 
M = 16. 



SNR in the observation channel is more pronounced compared 
to the SNR in the transmission channel. The case of high 
SNR = 20 dB and low SNR C = 10 dB is preferred by the 
estimator compared to the case of low SNRo = 10 dB and 
high SNR C = 20 dB. 

A set of box plots showing dependence of the SE and the 
ISE on the number of sensors distributed over the area A 
for M = 16 and M = 32 have been also generated. The 
results are similar to those for the case of M = 8 with the 
difference that the number of outliers as a function of the 
threshold decays faster to zero (see Figs. V\_ 12 T3j and 14) 
for illustration). 



function of varying threshold for three different realizations 
of SNRq and SNR C . Note that for M = 8 the effect of the 



V. Summary 

In this paper, a distributed ML estimation procedure for es- 
timating a parametric physical field is formulated. An iterative 
linearized EM solution is presented and numerically evaluated. 
The model of the network assumed (1) independent Gaussian 
sensor and transmission noise; (2) quantization of sensory data 
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Fig. 13. A box plot of the SE between the estimated and true location of 
the object displayed as a function of the number of sensors distributed over 
the area A. The number of quantization levels is set to M = 32. 
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Fig. 14. Dependence of the simulated ISE on the number of sensors 
distributed over the area A. The number of quantization levels is set to 
M = 32. 



prior to transmission; and (3) parametric function estimation 
at the FC. The stability of the EM algorithm has been 
evaluated for three different values of SNRo and SNRc- 
The results show that for a small number of quantization 
levels (quantization error is large) SNRo dominates SNRc 
in terms of its effect on the performance of the estimator. 
Also, when the sensor network is sparse, K = 10, 15, 20 
the EM algorithm produces a substantial number of outliers. 
Denser networks, K > 20, are more stable in terms of reliable 
parameter estimation. A similar analysis has been performed 
for M = 16 and M = 32. 

In the future, we plan to analyze the estimation abilities of 
the network at low SNR values and develop a Cramer-Rao 
bound on the estimated parameters. 



Appendix A 

This section provides details leading to the equation ([5]). 
Consider the z-th term under the sum in ([?]): 
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Note that the difference (ri — Gi) in the last integral can be 
rewritten as (n - Gf } + Gf } - G»). Then 
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Replacing the last integral with a difference of two Q- 



functions Q 



T-j — G 



(fc) 



and Q ' T -^- G 



(fe) 



we obtain: 



E* 



(r< - GO 



9Gi 

06> t 



x m ex p(-(^) 5Gi 



EE 



(fc)^2^ 




G,=G< 



(fc+i) 



References 

[1] J. Li, and G. AlRegib, "Distributed Estimation in Energy-Contrained 
Wireless Sensor Networks," IEEE Trans, on Signal Processing, vol. 57, 
no. 10, pp. 3746-3758, 2009. 

[2] D. Schabus, T. Zemen, and M. Pucher, "Distributed Field Estimation 
Algorithms in Vehicular Sensor Networks," The 73rd IEEE Conf. on 
Vehicular Technology, pp. 1-5, July 2011. 

[3] A. Ribeiro, G. B. Giannakis, "Bandwidth-Constrained Distributed Esti- 
mation for Wireless Sensor Networks - Part I: Gaussian Case," IEEE 
Trans, on Signal Processing, vol. 54, no. 3, pp. 1131-1143, 2006. 

[4] R. Niu and P. K. Varshney, "Target Location Estimation in Sensor 
Networks With Quantized Data," IEEE Trans, on Signal Processing, vol. 
54, no. 12, pp. 4519-4528, December 2006. 

[5] J. Y. Wu, Q. Z. Huang, and T. S. Lee, "Energy-Constrained Decentral- 



ized Best-Linear-Unbiased Estimation via Partial Sensor Noise Variance 
Knowledge," IEEE Signal Processing letters, vol. 15, no. 4, pp. 33-36, 
2008. 

[6] S. Cui, J. J. Xiao, A. J. Goldsmith, Z. Q. Luo, and H. V. Poor, "Estimation 
Diversity and Energy Efficiency in Distributed Sensing," IEEE Trans, on 
Signal Processing, vol. 55, no. 9, pp. 4683-4695, 2007. 

[7] Y. Wang, P. Ishwar, and V. Saligrama, "One-Bit Distributed Sensing and 
Coding for Field Estimation in Sensor Networks," IEEE Trans, on Signal 
Processing, vol. 56, no. 9, pp. 4433-4445, 2008. 

[8] R. D. Nowak, "EM Algorithms for Density Estimation and Clustering in 
Sensor Networks," IEEE Trans, on Signal Processing, vol. 51, no. 8, pp. 
2245-2253, 2003. 

[9] A. P. Dempster, N. M. Laird, D. B. Rubin, "Maximum Likelihood from 
Incomplete Data via the EM Algorithm," Journal of the Royal Stat. Soc. 
Series B, vol. 39, no. 1, pp. 1-38, 1977. 



